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Abstract 

We study a simple one-dimensional model of ballisticaly-controlled annihilation in 
which the two annihilating species are initially spatially separated. The time depen- 
dent properties of the annihilation front are exactly derived. It is shown that the 
front wanders in a brownian fashion around its average value. 
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1 Introduction 



The study of the kinetics of diffusion-controlled or ballisticaly-controlled an- 
nihilation provides a nice playground on which the tools of nonequilibrium 
statistical mechanics can be tested [1,2]. Particularly in low dimensions, the 
kinetics of these systems is governed by the fluctuations and any mean-field 
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like (or Boltzmann-like) approximation will not give a satisfactory description 
of the properties of the system. 

A lot of efforts have be expended in studying the diffusion-controlled case, both 
for the situation in which the reactants are uniformely distributed in space [3] 
and for the one in which they are initially spatially separated, leading to 
reaction-diffusion fronts [4,5]. In both cases, the long time behavior of physical 
observables is characterized by power law. The associated exponents have some 
universal properties [2]. 

More recently, the case of ballisticaly-controUed annihilation has been inves- 
tigated for spatially homogeneous conditions [6]. In particular, exact results 
have been obtained for one-dimensional, single species models by Piasecki [7] 
and Droz, Rey, Prachebourg and Piasecki [8,9]. 

In these models, one considers point particles which move freely, with a given 
velocity. When two particles collide they instantaneously annihilate each other 
and disappear from the system. The system with only two possible velocities 
-\-c or — c has been studied in a pioneering work by Elskens and Prisch [10]. 
The case of an arbitrary discrete velocity distributions has been investigated 
by Droz, Rey, Prachebourg and Piasecki [9]. An exact equation for the sur- 
vival probability until the time t of a particle moving with velocity v was 
derived, allowing to compute the density and the time dependent velocity dis- 
tribution in the asymptotic regime t ^ oo. Por a symmetric three velocities 
distribution, different kinetic regimes were found as a function of the initial 
fraction of particles at rest. Such processes can model several physical situa- 
tions as a recombination reaction in the gas phase or the fluorescence of laser 
excited gas atoms with quenching on contact (the one-dimensional aspect can 
be obtained by working in a suitable porous media [11]) or, the annihilation 
of kink-antikink pairs in solid state physics[12]. 

A related but different class of problems is the case in which the two species 
(called A and 5), are initially separated in space and annihilate on contact. 
Such process can model the situation in which chemical species incorporated 
in a gel move ballistically under the action of a drift [13]. As the two species 
cannot penetrate one into the other (as they annihilate on contact), a well 
defined reaction front is formed. We aim at computing exactly the time de- 
pendent properties of this front. 

The paper is organized as follows. In section two, the model is defined. In 
section 3, we compute the average number of collisions that have occured 
before time t. An exact analytic expression for the probability density to find 
the front at a given point and a given time is derived. In section 4, it is shown 
that the front makes a random walk around its (time dependent) average 
value. Possible extensions of this work are discussed in the conclusion. 
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2 The model 



We consider a one-dimensional system formed of two species of particles. Ini- 
tially, particles A are spatially randomly distributed in the region (—00, 0) 
and the B ones are spatially randomly distributed in the region (0, 00). The 
positions of the particles obey a Poissonian distribution. The initial linear 
average densities for the A and B particles are respectively and p^. The 
velocities of each particle is an independent random variable taking the value 
±c with even probabihty. Particles of the same kind suffer elastic collisions. 
When two particles A and B meet, they annihilate. Thus, practically the A 
particles with velocity — c and the B particles with velocity +c move freely. 
Accordingly, the relevant part of the dynamics concerns the A particles with 
velocity +c and the B ones with velocity — c. 

Let {yi,y2, ■ ■ ■ ,yk, ■ ■ ■) be the initial positions of the A particles with velocity 

+c and {xi,X2, ■ . . ,Xk, ■ ■ ■) the initial positions of the i?'s with velocity — c 
(see figure 1). The relative velocity between the A and B particles being 2c, 
the pair of particles initially at {xk,yk) will collide at time tk = {x^ — yk)/'^- 



A particles ^ B particles 

Vk ^ ^ X 2 Xk 

Vk-i ■■■ ys yi xi X3 ■■■ Xk^i 

Fig. 1. Typical initial condition. The particles are labeled such that 

Vk < ■ ■ ■ < yi < < xi < ■ ■ ■ < xk- 



The position at which this collision will take place defines the position of the 
annihilation front. Thus, the position of the front at time t is — yi), if 
2ct < Xi + yi (i.e. if no particle have yet collided), ^(x2 — 1/2), if xi + yi < 
2ct < X2 + y2 (i-e. only the first particles on the right and on the left have 
collided), and so on. The dynamics in itself is purely deterministic, the only 
stochastic aspect comes from the initial conditions. The properties of the front 
are completely defined by the probability density iJ,{X;t) to find the front at 
the point X, at time t. Moreover, another quantity of interest is N{t), the 
number of particles that have been annihilated until a given time t. 
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3 Analytical solution 



We first investigate the behavior of N{t). For this purpose, we introduce the 
probabihty P{k;t) that exactly k collisions have occured within the time in- 
terval [0,t]. For a given initial configuration, at least k collisions have occured 
before the time t, if 2ct > x^—yk, and at most k collisions have occured before 
t, if 2ct < Xfc+i — Hk+i, with Xk < Xk+i and yk > Vk+i- Hence averaging over 
the initial conditions, we find: 



P{k; t) = {9 {2ct - [xk -Vkl) 9 {[xk+i - yk+i] - 2ct)) 

= {9 {2ct - [xk - Vk])) - {9 {2ct - [xk+i - Vk+i])) ■ (1) 

where 6 is the usual Heaviside function and the brackets denote the average 
over the initial conditions (see Appendix A). Thus: 



N{t) = J2kP{k;t), (2) 

k=l 



which, using (1), reduces to 



N(t)^Y.(H'^ct-[xk-yk])). (3) 
fe=i 



Thus, using the results of Appendix A, N{t) reads: 



k=l ^ n V /■ 



R 

Noticing that the sum over k can be exchanged with the two integrals, we find 
(where a ^ ^{xk - yu) and h^\{xk + yu)) 

ct a 

N{t) = 2 p^p"" J da J dbe-^(^^+''^)e^(''^-'''')lo (2^pA^B - b2)) , 
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Iq is the modified Bcsscl function. Tlie integral over b is calculated in ap- 
pendix B and the integration on a is straightforward. We obtain 



2pV 



ct - 



2(pA + pB) 



(4) 



In terms of the dimensionless time 



T ^ Ct{p^ + p^) 



and the parameter 



A^B 



A = 



4 p^p 



we have: 



iV(r) = \ 



(8) 



After a transient regime, N{t) grows linearly with time. 

The front position is defined as the point where the next collision will take 
place. Thus the probability density p{X;t) to find the front at the point X, 
at time t is given by 



p{X; t) e{x^ - t/i - 2ct) 

+ /£ 5 fx - ^^±1±^) e {2ct - [x, -y,])9 {[x,+, + y^+i] - 2ct) 



\k=l 



6 X- 



xi +yi 



The first term in the right hand side reads simply 



(6) 

e{2ct-[xk-yk\)) ■ 



5[X 



xi + yi 



CO U 

= p^p^ j dxi j dyi e-"^"! e-^^yi 5 (^X 

— oo 

[e{X) e-2p^^+^(-X)e2^^^ 



xi +yi 







(7) 
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To calculate the second term in the right hand side, we remark that one has 

to average over the position x^., x^-^-i, yu and yk+i- We shall then first average 
over particles + 1, to obtain an expression containing only the position 
and yfe. To do this, we need to calculate the following integral 

Xk — oo 



putting a = U^k+i - Vk+i) and b = l{xk+i + Vk+i), we have: 



oo oo 



S^2p^p^ J da J dhe(a, + h-x^) e {a -b + yk)5{x-b) 



— oo — oo 



X e 



(pA+pB)„+(pA_^B)6 



The integral over b is straightforward, and the integral over a gives after some 
algebra 



A^B 



2pV 



Xk + yk 



Hence, the second term in the right hand side of eq. (6) becomes 



f xk + Vk 
\ 2 



S[X-'^^l±y^We{2ct-[xu-yk])) . 



Now using the relation (see appendix A) 



J2(f(^k + yk:Xk-yk)) 



k=l 



CO a 



(8) 



:2pV/ da I dbe-(''^+''^^^+(''^-''^)^f(2b,2a)Io (2^Jp^p^ (a? - b^)) , 



-a 

one finally finds: 
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fiiX;t) 



2pV 

ct a 



(9) 



2pV/ da J dbe-(''^+''^)^+(''^-''^)Mo (^2^pV (a^ - b^) 



A^B 



X -6{X -b) + 



p^ + p^ 



6) 



In terms of the dimensionless variables r and A, and the dimensionless coor- 
dinate ^ = {p^ + p^)X, we can write 



T a 



(10) 



with 



p^ + p^ 



Note that the function between the curly brackets is an even function of ^; we 
shall call it r). The rest of the paper will be devoted to the study of this 
function. 



4 Results and interpretation 



Starting from (10) and using the integrals quoted in appendix B, it is possible 
to verify that p{x; t) is correctly normalized to 1. Moreover, we can also calcu- 
late its different moments. However for this purpose, it is simpler to consider 
the Fourier transform of p. Let us define the function /i(p; r) as: 



p{p;t)= I dXe-(''^+'''')P^e-(''^-''^)^/i(X;r) 

— oo 

1 °° 

^-j dee-P«/i(e;r). (11) 

— oo 
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Here again, this integral can be calculated using the formulas quoted in ap- 
pendix B. One finds: 



fl{p;T) = Ae" 



1 



l-p2 



cosh (^^Jp^ + Xt^ + 



sinh (^y/p^ + Xt 



[12) 



It is now easy to calculate the different moment of n, using the fact that 



dX 



(p^ + p^)XXX;t) = (-ir 



■ (13) 



P= — K 



Noticing that k^ + A = 1, it is straightforward to check the normalization. The 
first moment ^(r) = (^) is given by 



f (r) = K 



2k 



[1 + N{r)] , 



(14) 



and the variance is 



4 2 3 5A , 
A^ A 4 4 



(15) 



A---2r(l-A) 



e-2--i^e-^- 



If p^j^p^, and in the long time limit, the front moves with a velocity pro- 
portionnal to the density difference. The case A = 1 is of particular interest, 
because it corresponds to the symetric case p^=p^. As expected: 



e(r) ^ 



(16) 



but 



4LW = 2[1 + A^(r)] 



(17) 



In the long time limit, we see that ■Crmsl''") ^he front move essentially as an 
unbiased random walker. Thus, we can expect p to be essentially a gaussian. 
In fact, we shall prove that for any A, the long time limit of is a gaussian. 
Starting from the definition of an arbitrary moment of p (eq. (13)), we may 
write for long time: 



(n = (-l)"Ae-- 



QU 

9p" 



(18) 



P= — K 
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using 



Qn 

dp"" 



1^ 



Qn 
Q^n 



+ 0(r 



(19) 



we obtain = {ktY + 0(r"-i). 

If A < 1, K 7^ and we know the leading term for every moment of If A = 1, 
the situation is different and will be considered later. 

It is known that for a gaussian distribution with a non vanishing mean ^(r) 
and a variance Crms (''')> 



gauss 



Srms V ' / I 



,t^o(n-2fc)!fc! \2e{T) 



(20) 



In the long time limit: 



£^ (r) 

Srms V / 

2e(r) - 



(21) 



and (e) 



gauss 



(KT)" + 0(r"-i). 



Thus, when A < 1 and for sufficiently long times (i.e. when the corrections 
to the leading term become negligible), the probabilty density takes the 
gaussian form: 



/^(C; t) « exp 



(e - ^r) 
2Ar 



(22) 



The front wanders around its mean value with an amplitude of order of mag- 
nitude \/Ar, like a biased random walker with a diffusion coefficient of A/2. 

The A = 1 case is different because the amplitude of the r" term vanishes. 
However, we can still prove that for large time, the probability density /i is a 
gaussian. More precisely, we shall give an upper and a lower bound, both of 
gaussian form. For this purpose, we shall come back to the Fourier transform 
ji{p\ r) defined by cq. (11). It can be shown (see appendix C), that the inverse 
Fourier transform of fl{p, r) is given by 



/i(^;T) = Ae 



, , X sinh ( vT+Ar 

cosh VI + Ar H \ — 

^ ' VI + A 



(23) 
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I . , . cosh (Vl + At - l^l) 

- \e{T - <j sinh (x/TTAt - lel) + ^^^^^ ^ 



- 2 



dr cos(r^) 
TT 1 + 



sinh {\J A — Hr^ + 



cosh ( VA — r^r 



\/A^ 



> . 



It is now easy to find an upper and a lower bound for /i. Indeed, for |^| < 
one has: 



cos(r^) ^ cos(r.^) 
1 + A - l + r2 



< cos(r^) 



as r e [0, \/A]. The remaining integral can be calculated exactly: 



cosh {\/\ — r^r 



— / drcos(r^) sinh ( vA — r^rj H ^ 

^ L ^ 

h (^/A(T2-e)) / / \ 



(24) 



yielding, then two bounds for /x. Two conclusions can then be drawn. First for 
A = 1, the asymptotic time behavior of (24) is 



TTT 



and then, for ^ < | and r large, we have: 



TTT 



TTT 



(25) 



In the long time limit, the density probability \x becomes the density proba- 
bility of a random walk, with a diffusion coefficient of |. 

The second consequence concerns the case ^ = 0, for any value of A: we have, 
in the long time limit: 



g(l-v/A)r 

/i(0; r) oc — — if A < 1 , 

/t(0;T)oc— p if A = 1 . 



(26) 
(27) 
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These results are not really surprising. Indeed, it can be verified that ^ = is 
a maximum of /i, for A = 1 but not for A < 1. 

It is useful to get an image of what is going on by plotting the function 
for various value of A. However for this purpose, both equations (10) and (23) 
are inadequate for technical reasons. Equation (10) is not well suited because 
it contains a double integral which requires a lot of CPU time and equation 
(23) because the integral diverge for r = \/A, leading to numerical difficulties. 
To avoid these problems, we have obtained a different integral representation, 
starting from (10), and using the following relation: 

h [^/X{r' - e)) = A/o (\/A(r2 - e)) ■ (28) 



This relation follows directly from the integral representation for Iq 



dz 1 
2i7r z 



exp 



2 



{a + P)z + 



a — P 



where C is, for example, the circle of radius 1 centered in z 
in (10), integrating twice by part and using that 



0. Inserting (28) 



= 2e-"5(/3-0, 



we get 



A(C;^) = A< 



e-l«l + lj da 



e-r |^e-l"-«l+e-l"+«l) 



Ji jA(r2-a2 

X I At ' +iJ^xiT^-a^) 



Performing the integral over the two theta functions, rearranging the terms in 
the remaining integral and using the integrals of appendix B, leads to: 
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Fig. 2. Plot of /i(C;r), for fA= 1 and p^= 1 (i.e. A = 1). 

/i(e;T)=Ae-l«l- 



cosh (Vl + At) + 



sinh ( Vl + At 



vm 



(29) 



A^(r-|e|)e-^ 



j da sinh(Q; — |^|) 



sinh(r — ^) 



1^1 



A(r2 - a2) 



+ Io 



A(r2 - 



We have not been able to simplify further this expression and in particular 

to perform analytically the integration. However, this form can be easily in- 
tegrated numerically. For examples, the function h{^]t) is plotted on figure 2 
and 3 for respectively A = 1 and for A = 0.5, p^>p^. 

One sees that for A = 1, the probability density spreads out and approches 
a gaussian for long times. For the assymetric case A = 0.5, one sees clearly the 
drift of the front which beecomes linear in the long time regime. 



5 Conclusions 



Using a combinatorial analysis approach, we have been able to derive exact 
results for the time dependent number of collisions, the probability density 
n{X;t) to find the front at point X at time t and its first two moments. 
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Fig. 3. Plot of /x(C;r), for p^= 1 + \^/?, and p^= 1 - \^/^ (i.e. A = 0.5). 

Although simple, this nonequilibrium model leads to nontrivial behavior as 
far as the dynamics of the annihilation front is concerned. 

Several extensions of the present work are possible. One is to consider a simi- 
lar model in which the velocities are not restricted to be bimodal, but obey a 
more general (discrete or continuous) distribution. Another interesting aspect 
is given by the generalization of this model to higher dimensions. In this case, 
some A particles can invade the right part of the system whithout being anni- 
hilated (and vice- versa for the 5's). A new definition of the front in term of the 
annihilation rate should be introduced and different time dependent behavior 
may occur. These more complicated problems are under investigation. 



A Appendix 



The average over the initial conditions is taken as follows: the B particles are 
distributed to the right of the origin, with a Poissonian distribution of density 
p^. Thus the probabihty to find the first particle at position xi is 

7ri(a;i,p^)=e-^'-ip^ . 
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The probability to find k — 1 particles between and Xk and a particle in Xk 
is 



. k-l 



The average of a function depending only on Xk is thus 

, k-l 



oo / I 

(fix,)) = J dxkf(xk) e-''"^^ Yk-i)! ■ ^^-^ 



Now, the average of a function depending both on Xk and Xk+i is 



oo oo 



{fixk,Xk+i))^ J dxk J dx,+if(xk,x,+i)e-"""^+^ ^a-i)! ^""^ " ^^"^^ 



Xk 



(k 



Indeed, the probability to find exactly k — 1 particles between and Xk, one 
in Xk and another in Xk+i, but none between is 



fc-1 



P e 



(^-1)! 



as e P^i^k+i ^k) pB ^YiQ probability that the interval [xk-,Xk+\] is empty 
Similar considerations apply for the A particles. 

Now, we propose to prove the formula (8). 



^{f{.Xk + yk,Xk-yk)) 



k=l 



7 r 

^Y^P^P^ J dxk J dykC 



Xk (-yk) 



k-l 



oo oo 



:2pV / da y db^(a + b)^(a-b)e 



[(k-l)!f 

xfi^Ck + Vk^Xk - Vk) 

(pA+pB)a+(pA_^B)b 



-OO — oo 



2pV/ da I dbe-(''^+''>+(''^-''^)'^f(2b,2a) 



oo a 



-a 



(A.3) 



(A.4) 
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where in (A. 4), we put 2b — Xk + Uk and 2a — Xk — yk and use the serie 
representation of Iq. 



B Appendix 



In this appendix, we quote two integrals that we often use in our derivation. 
The first is used, for example, in deriving the formula for N{t) (4), or in calcu- 
lating the different moments of without introducing the Fourier transform. 



y y 

- I dxe'^i- lo (/32v/y'-x^) = / dxcosh(/3ix)Io (^h^/7^) (B.l) 



2 

-y 



(see for example [14] equation 8). 

In calculating (28), we used the previous integral and also the following 



(B.2) 



Ii (yA(r2-a2)) 

/ dQ;cosh(Q;K)Ar — . = cosh(T) — cosh(ACT) , (B.3) 



/A(t2 - a2) 
which can be found in [14] equation 6) 



C Appendix 



In this appendix, we shall prove formula (23), starting from eq. (12). Remem- 
bering the definition (11) of /i(p;r), we may write 



/i(e;r) = Ae-- / ^e' 



p+ 1 p — 1 



X 



cosh ^y^p2 _|_ _|_ 



sinh (Vp^ + '^'7") 
VF+A 
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The integrant has two poles (p = ±1) and a cut on the imaginary axis between 
— i\/A and i^/X. For ,^ > r, we can close the contour to the left, without 
changing the value of the integral. We get then a contribution from the pole 
p — —1 and eventually from the cut (depending on the sign of k). Nevertheless, 
as cosh(x) +smh{x)/x is an even function, the contribution of the cut vanishes 
and we are left with 



/i(e;T) = Ae---« 



. . sinh ( Vl + Ar 

cosh I VI + At) H \ — 

^ ^ VI + A 



e>T. (C.l) 



For ^ < T, we close the contour to the right, leading to a contribution of the 
pole p — +1 only, as, again, the contribution of the cut vanishes. We find 



Me;T) = Ae— i«i 



cosh (Vl + Ar) + 



sinh ( vT+Ar 



VTTa 



e<-r.(C.2) 



For ^ e (— r, r), the situation is more complicated. First we write 



cosh ^y'p^ + Ar^ + 



1 + 



sinh + Ar^ 
Vp^TA 

1 \ -Vp^+At 



+ e- 



Vp^TA, 



For the first exponential, we can close the contour to the left and to the right 
for the second. Now, for both integrals, the integrant is not an even function 
of y/p'^ + A and one shall get a contribution of the cut. Integrating, one then 
eventually arrives to the formula 



/i(^;T) = Ae ^ < 



1 - 



v^r+A, 



cosh(^) 



(C.3) 



+ 2 



dr cos(r^) 
TT 1 + 



. . cosh ( \/A" 

sinh (^v A — r^r j H 



ie|<r. 



Combining (C.1,C.2,C.3) together gives (23). 
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